/***************************************************************************
* Piotr's Image&Video Toolbox      Version 2.2
* Copyright 2009 Piotr Dollar.  [pdollar-at-caltech.edu]
* Please email me if you find bugs, or have suggestions or questions!
* Licensed under the Lesser GPL [see external/lgpl.txt]
***************************************************************************/
#include "mex.h"

/***************************************************************************
* Return index of bin for x. The edges are determined by the (nBins+1)
* element vector edges. Returns an integer value k in [0,nBins-1]
* representing the bin x falls into, or k==nBins if x does not fall
* into any bin. if edges[k] <= x < edges[k+1], then x falls
* into bin k (k<nBins). Additionally, if x==edges[nBins], then x falls
* into bin k=nBins-1. Eventually, all values where k==nBins should be ingored.
* Adapted from \MATLAB6p5\toolbox\matlab\datafun\histc.c
***************************************************************************/
int findBin( double x, double *edges, int nBins )
{
	int k = nBins; /* NOBIN */
	int k0 = 0; int k1 = nBins;
	if( x >= edges[0] && x < edges[nBins] ) {
		k = (k0+k1)/2;
		while( k0 < k1-1 ) {
			if(x >= edges[k]) k0 = k;
			else k1 = k;
			k = (k0+k1)/2;
		}
		k = k0;
	}
	/* check for special case */
	if(x == edges[nBins]) k = nBins-1;
	return k;
}

void assignToBins( double *B, double* A, double* edges, int n, int nBins )
{
	int j; for( j=0; j<n; j++ ) B[j]=(double) findBin( A[j], edges, nBins );
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	int n, nBins; double *A, *B, *edges;

	/* Error checking on arguments */
	if( nrhs != 2 ) mexErrMsgTxt("Two input arguments required.");
	if( nlhs > 1 ) mexErrMsgTxt("Too many output arguments.");

	/* extract inputs */
	A = (double*) mxGetData(prhs[0]);
	n = mxGetNumberOfElements(prhs[0]);
	edges = (double*) mxGetData(prhs[1]);
	nBins = (int) mxGetNumberOfElements(prhs[1]) -1;

	/* create outputs */
	plhs[0] = mxCreateDoubleMatrix(1, n, mxREAL );
	B = (double*) mxGetData(plhs[0]);

	/* calculate the histograms */
	assignToBins( B, A, edges, n, nBins );
}
